Charmonium at high temperature in two-flavor QCD 
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We compute charmonium spectral functions in 2-fiavor QCD on anisotropic lattices using the 
maximum entropy method. Our results suggest that the S-waves {J/ip and 77c) survive up to tem- 
peratures close to 2T C , while the P-waves (xco and Xci) melt away below 1.2T C . 



I. INTRODUCTION 

The properties of hadrons or hadronic resonances 
above the deconfinement transition is a subject at the 
heart of the current experimental programme at RHIC. 
Questions of interest include the issue of which hadrons 
survive as bound states in the quark-gluon plasma, and 
up to which temperature; as well as the transport prop- 
erties of light and heavy quarks in the plasma. 

Of particular interest are charmonium states, follow- 
ing the suggestion Q that J/ip suppression could be a 
probe of deconfinement. Potential model calculations us- 
ing the heavy quark free energy have tended to support 
this picture. However, previous lattice simulations in the 
quenched approximation 0, [H, 0, H| indicate that con- 
trary to this, J /ib may survive up to temperatures as 
high as 1.5 — 2T C . Recently, potential model calculations 
using the internal energy of the heavy-quark pair have 
reached the same conclusion, and using the most recent 
lattice data [f| these models indicate a qualitatively sim- 
ilar picture in the case of N f =2 QCD [1,0,1. Support 
has also been provided by studies employing a real-time 
static potential U3] and a T-matrix approach which in- 
cludes scattering states [ll|. Note, however, that doubts 
have been expressed whether any potential model can 
satisfactorily describe the high-temperature behaviour of 
quarkonium correlators [12j , while some recent potential 
model studies have questioned the survival of quarkonia 

There are now high-statistics data available for J/tp 
production at SPS [llll5j and RHIC [Hj], showing simi- 
lar amounts of suppression at both experiments, despite 
the big difference in energy density. Two different sce- 
narios have been developed to explain this result. The 
sequential suppression scenario [l7| takes its cue from 
lattice results, suggesting that the entire observed sup- 
pression originates from feed-down from the excited IP 
and 2S states, which melt shortly above T c , while the IS 
state survives up to energy densities higher than those 
reached in current experiments. On the other hand, the 
regeneration scenario [H, [lj| HtJ HH suggests that ad- 
ditional charmonium is produced at RHIC energies from 



coalescence of c and c quarks originating from different 
pairs. 

Lattice simulations with dynamical fcrmions (2 or 2+1 
light flavors) will be one of the essential ingredients in 
resolving several of these issues. In the present paper, we 
present first results from such simulations. Preliminary 
results have appeared in Refs. [22l.[23|. 

Hadron properties are encoded in the spectral func- 
tions pr(io,p), which are related to the imaginary-time 
correlator Gy(r, p) according to 



Gr(r,p) 



2^ 



(1) 



where the subscript T corresponds to the different quan- 
tum numbers. The kernel K is given by 



K(t,w) 



cosh[w(r - 1/2T)] 
sinh(w/2T) 



(2) 
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From now on we consider zero momentum only and drop 
the p dependence. 

Spectral functions can be extracted from lattice corre- 
lators G(t) using the Maximum Entropy Method (MEM) 
(2^ . For this to work and give reliable results, it is neces- 
sary to have a sufficient number of points in the euclidean 
time direction: at least O(10) independent lattice points 
are needed. At T ~ 2T C , this implies a temporal lat- 
tice spacing a T < 0.025 fm. If the spatial lattice spacing 
a s were to be the same, a simulation with dynamical 
fermions on a reasonable volume would be far too expen- 
sive to carry out with current computing resources. 

In order to make the simulation feasible, anisotropic 
lattices, with a T <C a s , are therefore required. However, 
dynamical anisotropic lattice simulations introduce addi- 
tional complications not present in isotropic or quenched 
anisotropic simulations. The anisotropic formulation in- 
troduces two additional parameters, the bare quark and 
gluon anisotropies, which must be tuned so that the 
physical anisotropies are the same for gauge and fermion 
fields. In the presence of dynamical fermions, this re- 
quires a simultaneous two-dimensional tuning, which has 
been described and carried out in Ref. [25| . 

In this study we attempt to determine charmonium 
spectral functions in 2-flavor QCD using anisotropic lat- 
tices and the Maximum Entropy Method. The MEM 
analysis has been performed using Bryan's algorithm 12611 
with the modified kernel recently introduced in Ref. [27| . 
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We found that this greatly improved the stability and 
convergence properties of MEM. In Sec. |TT] we describe 
our procedure and simulation parameters. In Sec. IIIII we 
briefly discuss the spectrum at zero temperature, while 
Sec. lIVI contains the main body of our results above T c . A 
detailed discussion of dependence on the default model, 
time range, energy cutoff and statistics is given in Sec.lVl 
Finally, in Sec. I VII we discuss remaining uncertainties and 
give our conclusions and prospects for further work. 

II. SIMULATION DETAILS 

We use the Two-plaquette Symanzik Improved gauge 
action [28| and the fine-Wilson, coarse-Hamber-Wu 
fermion action [29^ with stout-link smearing [3Cj |. The 
process of tuning the action parameters, and the param- 
eters used, are described in more detail in Ref. [25|. We 
have performed simulations with parameters correspond- 
ing to Run 6 in Ref. [25J as well as at the tuned point, 
which we denote Run 7. The parameters are given in 
Tables U and [TTJ They correspond to a spatial lattice 
spacing a s w 0.165 fm with a (renormalised) anisotropy 
£ = a s /a T w 6. The sea quark mass corresponds to 
rrin/mp w 0.54. The lattice spacing was determined from 
the IP-IS splitting on the 12 3 x 80 Run 7 lattice; the Run 
6 lattice spacing was then determined using the IP-IS 
splitting on the 8 3 x 80 lattice corrected for finite volume 
effects. 

The pseudocritical temperature T c was determined by 
measuring the Polyakov loop (TrL) on 12 3 x N T lat- 
tices on Run 6. A jump in the value of (TrL) was 
found between iV r = 34 and 33, so we conclude that 
a T T c « 1/33.5, or 205-210 MeV for both parameter sets. 
We have not been able to determine the pseudocritical 
temperature T c to greater precision on these lattices be- 
cause of the finite lattice size. Partly for this reason, we 
have chosen to express our temperatures in units of MeV 
rather than as T/T c , as is often done in quenched stud- 
ies. Since this analysis is carried out with 2 dynamical 
light quark flavors, there is also less need to rescale tem- 
peratures with T c to correct for the difference between 
the simulation and the real world with 2 + 1 light quark 
flavors. 

We have computed charmonium correlators in the 
pseudoscalar (77c) and vector (J/ip) channels, as well as 
the scalar (Xco) and axial-vector (xci) channels. In the 
nonrelativistic quark model, the former two are S-waves 
and the latter two P-waves. In this study we have used 
local (unsmeared) operators, 

Gr(r) = ww E< M ^ x '*) Mr (y' t+r ))' ( 3 ) 

s T x,y,t 

where 

M r (x,r) =^(x,r)r^(x,r), (4) 

and r = 75, 7^, 1, 757; for the pseudoscalar, vector, scalar 
and axial- vector channel respectively. All-to-all propaga- 



Run 






£ 9 £s a T 1 a s 


£° a T rn° c 


6 


8.06 


7.52 


5.90 6.21 7.06GeV 0.167fm 


5.9 0.08, 0.092 


7 


8.42 


7.43 


6.04 5.84 7.23GeV 0.163fm 


5.9 0.117 



TABLE I: Simulation parameters. c are the bare (input) 
anisotropics for gluons (g), sea quarks (s) and charm quarks 
(c), while £ 9l s are the renormalised (measured) anisotropies. 
The charm quark anisotropy was tuned independently to give 
an output anisotropy of 6. a T and a s are the temporal and 
spatial lattice spacings. The bare sea quark mass is a T m s — 
—0.057 for both sets of parameters, with m^/m p = 0.54. 



Run N B N T T (MeV) T/T c N cig 


6 8 


80 


88 


0.42 


250 


12 


33 


214 


1.02 


80 


8 


32 


221 


1.05 


500 


12 


32 


221 


1.05 


400 


12 


31 


228 


1.08 


100 


12 


30 


235 


1.12 


100 


12 


29 


243 


1.16 


100 


12 


28 


252 


1.20 


125 


8 


24 


294 


1.40 


1000 


12 


24 


294 


1.40 


500 


8 


20 


353 


1.68 


1000 


12 


20 


353 


1.68 


1000 


8 


18 


392 


1.86 


1000 


8 


16 


441 


2.09 


1000 


12 


16 


441 


2.09 


500 


7 8 


80 


90 


0.42 


250 


12 


80 


90 


0.42 


250 


8 


32 


226 


1.05 


1000 


8 


24 


301 


1.40 


1000 


8 


16 


451 


2.09 


1000 



TABLE II: Lattices and parameters used in this simulation. 
The separation between configurations is 10 HMC trajecto- 
ries, except for the N T = 80 runs where configurations were 
separated by 5 trajectories. 



tors [3l| have been used to improve the signal and sample 
information from the entire lattice. The propagators were 
constructed with no eigenvectors and two noise vectors 
diluted in time, color and even/odd in space. On the 8 3 
lattices, for Run 6, we have computed correlators for two 
different bare quark masses, a T m c — 0.080 and 0.092, as 
the precise charm quark mass had not been determined 
on these lattices. Both masses are somewhat smaller than 
the physical charm quark mass. This may result in an un- 
derestimate of the melting temperatures for the P-waves. 
For Run 7 we used a T m c — 0.117, tuned to reproduce the 
physical J/ip mass on the 12 3 x 80 lattices. The bilin- 
ear operators have not been renormalised, so our results 
only concern the shapes of the resulting correlators and 
spectral functions, not their overall magnitude. 
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Run 


a T m c 


mps 


my 


mAV rase 


6 


0.080 


2.643 


2.689 


3.118 3.018 




0.092 


2.800 


2.835 


3.233 3.209 


7 


0.117 


3.145 


3.174 


3.637 3.615 



TABLE III: Ground state masses (in GeV) at zero tempera- 
ture from a variational calculation. The a T m c = 0.08 results 
were obtained by extrapolation from two higher masses. 
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FIG. 1: Pseudoscalar spectral function at zero temperature 
on the 8 3 x 80 lattice (Run 7). The dashed line denotes the 
standard spectroscopy result quoted in Table ITTT1 



III. ZERO TEMPERATURE 

The charmonium spectrum at zero temperature (N T — 
80) has been computed using standard spectroscopic 
methods, with a variational basis of smeared operators 
in S-, P- and D-wave channels. Preliminary results were 
presented in Refs. l32l . [33|; the full results will be re- 
ported elsewhere |34| . Here we only report results for 
ground state S-wave (pseudoscalar, vector) and P-wave 
(axial, scalar) masses, which are given in Tabic [TTTT 

In Figure Q] we show the pseudoscalar spectral function 
for our T = lattice (8 3 x 80, Run 7). Each spectral 
feature is fitted to a Gaussian with peak position M, full 
width at half maximum T. The "error" bars shown in the 
figure require careful interpretation. The horizontal bar's 
centre and width represent M and T respectively, and its 
height represents the area of the Gaussian evaluated over 
the range M - T/2 to M + T/2. The vertical error bar 
represents the error in this area as determined by the 
Bryan algorithm (26j. The width of the horizontal bar 
does not correspond to the error in the peak's position. 
We expect that this width is primarily determined by 
statistics, and will decrease as our correlators become 
better determined, see Sec. [V] 

The position of the primary peak can be seen to agree 
with the standard spectroscopy results within errors. The 
second peak at 4.1 GeV cannot be identified with the first 
radial excitation r/ c (2S), which has a mass of 3.64 GeV; 
rather, it is most likely a combination of the 2S, 3S and 



4S states, possibly with some contamination from lattice 
artefacts. With more statistics it should be possible to 
resolve these states further, as has been demonstrated in 
quenched QCD some time ago [35|. The third bump in 
the spectral function is most likely a lattice artefact, cor- 
responding to a cusp in the free lattice spectral function. 
As shown in the Appendix, the free spectral function has 
cusps at a T uj ~ 0.72 and 1.14, corresponding to 5 and 
8 GeV respectively; these may merge or be pushed to 
higher energies in the interacting case. 

We find the same picture in the vector channel. In the 
axial and scalar channels the spectral function is much 
less well determined; however, the position of the primary 
peak is found to agree within errors with the standard 
spectroscopy result also in these channels. 



IV. HIGH TEMPERATURE 

Spectral functions just above T c (T = 226 MeV, 
T/T c = 1.05) are presented in Fig. [2] We show results 
in four channels, on the 8 3 x 32 lattice (Run 7). To ob- 
tain these results, we used the continuum free spectral 
function m(u)) — m^uj 2 as default model and discretised 
the energy integral JT]) using a T Aw = 0.005 and a cutoff 
5.0 (w m ax = 35 GeV). Since the first two times- 
lices may contain short-distance lattice artefacts we have 
used G(t) at r/a T = 2, . . . , N T /2 in Eq. |T]). An extensive 
discussion on the dependence on these choices is given in 
Sec. |V] In all channels we find a peak which is consis- 
tent with the zero-temperature ground state mass. There 
are indications that the vector, axial-vector and scalar 
masses have shifted slightly upwards, although this can- 
not be determined with any certainty given our current 
precision. The second peak at u ~ 6 GeV is again most 
likely a lattice artefact, as discussed in the Appendix for 
the free theory. It should be noted that the proximity of 
this second peak may distort the shape of the primary 




co (GeV) 

FIG. 2: Spectral functions on the 8 3 x 32 lattice (Run 7), 
in the pseudoscalar (PS), vector (V), axial- vector (AV) and 
scalar (SC) channels. 



FIG. 3: Reconstructed correlator in the vector (J/VO and 
pseudoscalar (r/c) channel, for different temperatures, on 8 3 x 
N T lattices. The filled symbols are for Run 7, while the open 
symbols are for Run 6, a T m c = 0.092. 



FIG. 4: Reconstructed correlator in the scalar (Xco) and axial- 
vector (xd) channel, for different temperatures, on the 12 3 x 
N T lattice (Run 6, a T m c = 0.080). 



peak. In order to fully disentangle the first peak from 
any lattice distortions, simulations with finer lattices are 
necessary. However, at this temperature the structure in 
the spectral functions is quite robust and, given the po- 
sition of the first peaks, we are confident that they are 
separate features corresponding to the ground states in 
the respective channels. 



A. Reconstructed correlators 

One way of determining whether or not medium modi- 
fications of hadron properties are present, is by studying 
reconstructed correlators [3||. These are correlators ob- 
tained by integrating up Eq. |T]) with the spectral func- 
tion p(u>; To) obtained at some reference temperature To, 
and the temperature-dependent kernel K(t,lu;T). If the 
spectral function is unchanged, the reconstructed corre- 
lator GrecfV) will be equal to the actual correlator G(r), 
while, conversely, if G rcc (r) ^ G(t) this shows that the 
spectral function must be modified. In this procedure 
MEM is only used at the lowest temperature To (with 
the largest value of N T ), making this analysis a robust 
tool for higher temperatures. As we will demonstrate 
shortly, we find that the conclusions drawn from the re- 
constructed correlators in our dynamical simulations are 
surprisingly close to those obtained in quenched lattice 
QCD studies 0,i]. 

Figure [3] shows the reconstructed correlator in the S- 
wave (vector and pseudoscalar) channels, using the spec- 
tral functions obtained at T = 221 (Run 6) and 226 (Run 
7) MeV (N T = 32, see Fig. [5]) as the reference point. In 
the pseudoscalar channel we see very little change: only 
at the highest temperature (T = 441 and 451 MeV for 
Run 6 and 7 respectively; T/T c = 2.1) does the recon- 
structed correlator differ from the actual one by more 
than 3% at large imaginary time. This suggests that r\ c 



survives relatively unscathed in the medium up to this 
temperature, although it cannot be ruled out that even 
a 2% change in the reconstructed correlator may corre- 
spond to substantial modifications in the spectral func- 
tion [l3| . In the vector channel, somewhat larger modifi- 
cations are seen, although still less than 10% at the high- 
est temperatures. This may be related to the transport 
contribution which can be present in vector correlators, 
and is related to quark diffusion [53, [38, 39] . We have also 
compared the pseudoscalar correlator at N T — 32 with 
the reconstructed correlator from the zero-temperature 
spectral function shown in Fig. [TJ In that case we found 
no more than a 1.5% difference at large r. 

Figure 0] shows the reconstructed correlator in the 
P-wave (scalar and axial- vector) channels, again using 
T = 221 MeV as reference temperature. Here we see 
much greater changes in a smaller temperature range: al- 
ready at T = 235 MeV (T/T c = 1.12) the long-distance 
correlator differs from the reconstructed one by 20%, 
while at T = 252 MeV (T/T c = 1.2) the difference is 
up to 50%. If we instead use T = as reference temper- 
ature, we find that the T = 221 MeV correlator differs 
from the reconstructed one by a factor 2.5 at large dis- 
tances and by 20% at r/a T = 10. We infer that there 
are considerable medium modifications in this channel 
for T c < T < 1.2T C . Whether this corresponds to ther- 
mal broadening, a mass shift or melting of the \ci state, 
will be investigated in the following. 



B. Temperature-dependent spectral functions 

We now proceed to a discussion of temperature depen- 
dence of spectral functions in the range T c < T < 2.1T C . 
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FIG. 5: Pseudoscalar spectral function for different temper- 
atures on the 8 3 x N T lattice, for Run 7 (top) and Run 6, 
a T m c = 0.092 (bottom). All results have been obtained using 
m(uj) — 3lu 2 , tJmax = 35 GeV, and r/a r = 2, . . . , N T /2. 



1. Pseudoscalar channel 
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FIG. 6: As in the lower panel of Fig. [5] for the vector spectral 
function, using m(uj) = 8u> 2 . 



second peak, our spectral functions are nonzero every- 
where, and we are therefore not able to unambiguously 
distinguish the two possibilities. However, a threshold 
enhancement would be expected to become smaller as 
the temperature is increased, while we find a remarkably 
constant peak, consistent with a bound state. Because of 
these uncertainties, we are not in a position to conclude 
exactly when the rj c melts. However, our results suggest 
that the r\ c state is bound up to T w 392 MeV. 

In general, we see very little volume dependence in this 
channel, with the N s — 12 data for the most part being 
completely compatible with the N s — 8 data. This is 
consistent with r\ c being a compact bound state with a 
diameter much smaller than our lattice size, and indicates 
that this remains the case in the plasma up to the point 
where it melts. 



Figure [5] shows the pseudoscalar spectral function at 
various temperatures on the 8 3 x N T lattices. The r\ c 
peak can be seen to persist up to at least T = 392 MeV 
(A^ T = 18). At our highest temperature, T = 440 - 
450 MeV (N T = 16), no peak survives for Run 7 or for 
the larger lattice on Run 6, while the smaller lattice on 
Run 6 shows a distorted peak structure with a very large 
uncertainty in the peak height. Since the correlators on 
the two volumes differ by less than 0.5%, this discrepancy 
is more a sign of a breakdown of MEM than a physical 
effect. At these high temperatures the small number of 
available points means that it can not be determined at 
present whether the disappearance of the peak signals the 
melting of the resonance or the failure of the maximum 
entropy method. Indeed, the spectral function obtained 
from Run 6 N T = 32 correlators using the same time 
range (r = 2 — 8) and default model also exhibits no 
peak. 

The possibility that at higher temperatures there is no 
bound state, but only a threshold enhancement, must 
also be considered. Because of the proximity of the 



2. Vector channel 

The spectral function in the vector channel is shown 
in Fig. [HI We observe the same pattern as in the pseu- 
doscalar channel. The ground state peak appears to melt 
around 350 MeV (T/T c « 1.7, N T = 20), although it is 
again difficult to draw firmer conclusions, especially at 
higher temperature. At the highest temperature no peak 
is visible any more. Instead, we find nonzero spectral 
weight at all energies. This may be related to a trans- 
port contribution, signalling a nonzero charm diffusion 
coefficient. We hope to address this in the near future. 

3. Axial channel 

Fig. [7] shows the temperature dependence of the axial- 
vector spectral function on the 12 x N T lattice (Run 6, 
a T m c = 0.08). Since the P-waves are much more sensi- 
tive to finite volume effects than the S-waves, we use the 
larger volume in this analysis. The ground state peak 
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FIG. 7: Axial-vector spectral function for different temper- 
atures on the 12 3 x iV T lattice (Run 6, a T m c = 0.080). All 



results have been obtained using m(u) — 2ui 
GeV, and r/a T = l,...,N T /2. 



= 35 
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FIG. 8: As in Fig. [7] for the scalar spectral function, using 
m(uj) — u 2 . 



appears to survive up to T = 243 MeV (T/T c = 1.16, 
N T = 29), while at T = 252 MeV (T/T c = 1.2, N T = 28) 
there is no sign of any Xci peak. We interpret this as 
a sign of the melting of Xci somewhere in this tempera- 
ture range. A more detailed study of the 12 3 x 28 data 
reveals that by varying m(uj) or w max it is, however, pos- 
sible to reproduce a weak \ci peak, indicating that the 
bound state may still not have completely disappeared 
at this point. Higher statistics and lattices closer to the 
continuum limit will be required to resolve this issue. 



4- Scalar channel 

Finally, in Fig. [5] the scalar spectral function is shown 
for temperatures ranging from 219 MeV to 252 MeV. We 
see a similar pattern as in the axial channel, although 
the XcO state appears to melt at somewhat lower temper- 
ature than the Xci state: at T = 235 MeV (T/T c = 1.12, 



N T = 30) there is no sign of any surviving bound state. 
However, the scalar correlators are considerably noisier 
than the axial-vector correlators, so it is possible that we 
simply do not at present have sufficient statistics to ob- 
tain a signal in this channel. Indeed, given the slightly 
smaller change in the correlators shown in Fig.[4l a lower 
melting temperature seems surprising. Increased statis- 
tics will be required to resolve this issue. 



V. MEM SYSTEM ATICS 

In order to study the robustness of the spectral func- 
tions shown in the previous section, we now consider the 
dependence of the MEM output on the parameters that 
can be varied. This includes the default model depen- 
dence, dependence on the energy cutoff and discretisa- 
tion, dependence on the time range used in the analysis 
and the role of finite statistics. We focus on the pseu- 
doscalar and axial-vector spectral functions on lattices 
with a time extent of A^ r = 32, since we find that the 
vector and scalar channels behave qualitatively similar 
to the pseudoscalar and axial channels respectively 

We start with a discussion of the choice of default 
model. Since we are primarily interested in the prop- 
erties of the spectral functions in the 3 — 5 GeV region, 
we have mostly used the continuum free spectral function 
m(uj) = rriQUj 2 as default model, rather than the default 
model m(uj) = moOj(b + ui) proposed in Ref. (27j . which 
allows for nontrivial behaviour in the to — > limit. At 
the intermediate energies considered here, we find that 
the two models result in the same spectral function if 
the same value for the model parameter mo is used, al- 
though the second one tends to yield lower values for mo, 
when mo is determined by a single parameter fit to the 
correlator, using Eq. |T]). In addition, we have also used 
two other default models: m(u>) = mo and m(uj) = mow, 
which have very different high-energy behaviour. To as- 
sess the sensitivity of our results to the choice of default 
model, we have varied the parameter mo over a wide 
range. The output using these different models gives an 
indication of how tightly the data constrain the spectral 
function. 

Fig. [9] (top) shows the pseudoscalar spectral function 
for a large class of default models. The first three de- 
fault models vary in their normalisation over more than 
two orders of magnitude. Since the vertical axis of Fig. [5] 
is p(uj)/ui 2 , these three default models could be plotted 
as horizontal lines, at 0.3, 8, and 80 respectively. The 
fourth and the fifth default model differ from the first 
three particularly at small u>. The final two default mod- 
els {m{uj) — mo and m(uj) = mQLu) behave in a qualita- 
tively different manner, as 1/ui 2 and 1/lo respectively in 
this plot. In the absence of any input information from 
the euclidean correlators, the MEM output reproduces 
the default model. Since this is not happening here, 
we conclude that the MEM procedure is fairly robust 
against variations in the default model. In particular, 
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FIG. 9: Pseudoscalar spectral functions on the 8 3 x 32 lattice 
(Run 7) , for different default models (top) and energy cutoffs 
(bottom). 



— m(co) 


= 0.1o) 2 


— m(co) 


= 2.0o) 2 


-— m(co) 


= 8.00)2 


-- m(0)) 


= 0.5(G)+a) 2 ) 


— • m(0)) 


= 0.00066 


- m(co) 


= 0.1776a) 




4 6 

co (GeV) 



FIG. 10: Axial- vector spectral function on the 12 x 32 lattice 
(Run 6) with a T m c — 0.092, for different default models (top) 
and different energy cutoffs (bottom). 



we find that the leading edge of the spectral function is 
very robust, while also the height and position of the first 
peak are reasonably independent of the choice of default 
model. 

For some choices of default models parameters (espe- 
cially for larger values of mo) there appears to be a middle 
peak just above 4 GeV, or a broadening of the primary 
peak. This second peak, when it appears, coincides more 
or less with the second peak in the zero-temperature spec- 
tral function. This is too high to correspond directly to 
the radial excitation, rf c (3638MeV), but it might corre- 
spond to a radial excitation modified by medium effects 
and the nearby lattice doubler. However, since this peak 
is not reproduced for most of the parameters shown, we 
are cautious attaching too much physical value to it. 

The energy integral (JTJ has been discretised with 
a T Aw = 0.005 and a cutoff at w max . We have studied the 
sensitivity of the results to the cutoff by varying w max , 
while keeping Aw fixed; in practice we find that varying 
Auj does not change the results. In Fig. [9] (bottom) we 
show the dependence of the pseudoscalar spectral func- 
tion on the energy cutoff w max . We find little sensitivity, 
provided that w max > 28 GeV, or a T ui max > 4. 

We have performed the same analysis also on the Run 6 



lattices, for both charm quark masses and both volumes, 
and find very little dependence on either energy cutoff or 
default model in this case. 

In the axial-vector and scalar channel we expect finite 
volume effects to be significant. Therefore we will here 
analyse the larger lattice, 12 3 x 32. 

In the top panel of Fig. [TO] we show the effect of vary- 
ing the default model m(u>) on the axial-vector spectral 
functions. There is a great deal of variation, but in all 
cases we find either a ground state peak corresponding 
to the Xd state and a second peak at 6 — 7 GeV, or a 
broad structure encompassing the two, with a plateau in 
the middle. In this case, we cannot say with any con- 
fidence whether what we see is a bound state peak or 
a continuum threshold, but the presence of a structure 
near the \ c mass may indicate that Xd survives at this 
temperature, close to but just above T c , albeit possibly 
in a modified form. Generically, we find that the spec- 
tral function analysis is less robust for P-waves than for 
S-waves, which may be due to the local operators used 
in this study. 

In the lower panel of Fig. [TO] we show the effect of 
varying the energy cutoff w max on the axial- vector corre- 
lator. We see very little dependence on the cutoff, in the 
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FIG. 11: Pseudoscalar spectral function on the 8 x 32 lattice 
(Run 7) for different time ranges used in the MEM analysis: 
fixed T m in (top) and fixed r max (bottom). 



FIG. 12: Pseudoscalar spectral function at higher tempera- 
ture on the 8 3 x 24 (top) and 8 3 x 18 (bottom) lattice (Run 
6), for different default models. 



range shown here, but for lower energy cutoffs, w max J$ 28 
GeV, the peaks become more "washed out" . We take this 
as evidence that although the maximum energy for free 
fcrmions is a T uj max = 1.48, in the interacting theory the 
spectral function reaches higher energies, which must be 
included in the integral. 

The effect of varying the time range (r m i n , r max ) used in 
the MEM analysis is shown in Fig.lTTIfor the pseudoscalar 
correlators. We find a reasonable stability in our results 
as long as at least 10 data points are included; for T m ; n = 
2 or 3 even fewer points are required to reproduce the 
spectral function. 

We have carried out the same analysis at all tempera- 
tures, in order to try to clarify whether the presence or 
absence of a ground state peak is a physical effect or an 
artefact of the MEM. This is illustrated in Fig. [TH for 
the pseudoscalar channel at T = 294 MeV (8 3 x 24) and 
T = 392 MeV (8 3 x 18). We see evidence of a surviving 
ground state (j] c ) peak, but there is a quite strong de- 
pendence on both default model and energy cutoff, which 
becomes stronger as the temperature is increased. This 
means that our data are not sufficient to unambiguously 
determine whether the bound state survives at these tem- 
peratures, much less to say anything quantitative about 



changes to the spectral function. 

Finally, spectral functions reconstructed using MEM 
on a finite sample will always display a finite peak width, 
so the width of the peaks found here cannot be di- 
rectly interpreted as a thermal width of the correspond- 
ing mesonic states. (A further limitation is given by the 
finite resolution offered by the singular value decomposi- 
tion procedure used in our MEM analysis, but we believe 
we are not yet in this regime when N T = 32.) One may 
attempt to disentangle the unphysical statistical width 
from a possible physical thermal width by varying the 
number of configurations used. Specifically, if the shape 
of the spectral function is found to be independent of the 
statistics above a certain number of configurations, one 
can be more confident in the physical relevance of the 
results. 

In Fig. Q2] we show the pseudoscalar spectral function 
on the 8 3 x 32 lattice, obtained using different numbers of 
configurations. We see that as the number of configura- 
tions is reduced, the primary peak at first gets narrower, 
then remains approximately constant before broadening 
somewhat for the lowest statistics. The rather surprising 
initial narrowing may be related to the disappearance of 
the weak secondary peak discussed above, in which case 
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FIG. 13: Pseudoscalar spectral lunction on the 8 3 x 32 lat- 
tice for varying statistics, with u; max = 35 GeV, r/a T = 
2, . . . , N T /2 and m(uj) = 3a; 2 (top), m{uj) — l&ui 2 (bottom). 

it may be argued that the peak width for intermediate 
statistics is in fact a real thermal width. To test this hy- 
pothesis, we also show, in the bottom panel of Fig. Q2J 
the spectral functions obtained from the same data but 
with m(uj) = 16a-> 2 , where we already have seen that 
three peaks are produced. We see that as the statistics 
are reduced, the middle peak vanishes, but the primary 
peak remains unchanged. Only for very low statistics 
did we find a broadening. This lends some support to 
the hypothesis that the middle peak is indeed related to 
a surviving rj' c state. However, this result must be treated 
with caution because of the proximity of the lattice arte- 
fact peak at u> ~ 6 GeV. 

VI. DISCUSSION AND CONCLUSIONS 

We have computed charmonium correlators at a range 
of different temperatures on anisotropic lattices with two 
light sea quark flavors. We find that the S-wave (vector 
and pseudoscalar) correlators remain largely unchanged 
as the temperature is increased up to about twice the 
pseudocritical temperature, or 400 MeV. The P-wave cor- 
relators, on the other hand, exhibit substantial modifica- 



tions already between 220 and 250 MeV. This behaviour 
of the correlators is in good agreement with what has 
been found in quenched QCD studies [3, Hi S @- Us- 
ing the maximum entropy method to obtain the cor- 
responding spectral functions, our results indicate that 
the ground state S-wave peak survives largely unchanged 
up to T ~ 390 MeV, while at our highest temperature, 
T « 440 MeV, uncertainties in the MEM procedure pre- 
vents us from drawing any conclusion about the presence 
or absence of a ground state. In the axial- vector (P-wave) 
channel, we find that the ground state appears to melt 
between 240 and 250 MeV, although higher statistics will 
be needed to draw definite conclusions. The scalar me- 
son Xco appears to melt earlier, although this may be a 
function of limited statistics. Generically, we find that 
the spectral function analysis for S-waves is more robust 
than for P-waves, which may be related to the local oper- 
ators used to represent the meson states. There is some 
indication that a radial S-wave excitation may survive in 
the plasma phase just above T c , but it is premature to 
draw any conclusions about this. Again these results are 
in qualitative agreement with most corresponding calcu- 
lations in the quenched approximation [1,11, 0, [J. 

Our results appear to be compatible with the sequen- 
tial charmonium suppression scenario [l7T ] , which requires 
that S-waves melt at T > 2T C while P-waves melt close 
to T c . A simple hydrodynamical model calculation based 
on the sequential suppression picture [4(| gave melting 
temperatures of 2.1T C for the S-waves and 1.34T C for the 
P-waves and radial excitation. The former is quite com- 
patible with our results, while the latter appears quite 
high; however, given the simplicity of the model calcula- 
tion and the systematic uncertainties in our calculation, 
the discrepancy is relatively minor. 

There are several features of this calculation which 
must be improved before any firm, quantitative conclu- 
sion can be reached. The most important of these re- 
late to the use of a single, fairly coarse lattice spacing, 
a s 0.17 fm and a T « 0.028 fm. As a result, we are un- 
able to reach temperatures much beyond 2T C or 440 MeV, 
and our results at the highest temperatures are subject 
to uncertainty due to the small number of points in the 
imaginary-time direction. Furthermore, lattice artefacts 
at larger energies expected from free fermion calculations 
are close to the first peak representing the groundstate 
at lower temperatures, which complicates a straightfor- 
ward interpretation. A finer lattice would help overcome 
both of these problems. Simulations on finer lattices, 
bringing the main systematic uncertainties in this study 
under control, are currently underway. 

We also note that the fairly heavy sea quarks bring T c 
up from its physical value of 150-200 MeV [IE [H, EH , 
as does the absence of a third active flavor. Lighter sea 
quarks will also facilitate charmonium dissociation and 
thus bring down the melting temperature. Simulations 
with lighter sea quark masses are planned. 

In terms of addressing the experimental situation, two 
further developments are possible. Firstly, the RHIC 



10 



experiment corresponds to a small but nonzero baryon 
chemical potential, while the calculations presented here 
have been carried out at zero chemical potential. It 
would be useful to calculate the response of the meson 
correlators to a small chemical potential to determine 
what, if any, effect this has. Secondly, and perhaps more 
importantly, the J/ip particles which escape from the 
plasma and are observed as dileptons in detectors will 
have nonzero (transverse) momentum, and the momen- 
tum and rapidity dependence of the J/ip yields is a cru- 
cial factor in differentiating different models [l^, [2(| • It 
is therefore important to study the temperature depen- 
dence of charmonium correlators and spectral functions 
also at nonzero momentum. This is currently underway. 
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APPENDIX A: FREE LATTICE SPECTRAL 
FUNCTIONS 

In order to estimate lattice artefacts, we have studied 
meson spectral functions in the free lattice theory, follow- 



ing the approach of Refs. [H, H^, |4(| . Since the temporal 
discretisation in the fermion action used in this paper is 
identical to the standard Wilson fermions, most details 
can be found in Section 3.1 of Ref. (45|. Here we only list 
the expressions that are different. 

The fermion dispersion relation is determined by 

cosh(^ = l + ^J, (Al) 
where in this case 

3 

£k = 777 V* li (8 sin h - sin 2ki) , 
2s 3 

■Mk = ^ r a T Tn + — (3 — 4 cos ki + cos 2ki) , (A2) 
^ i=i 

with \i r = 1 + a T m/2 and s = 1/8. The free meson 
spectral functions take the same form as in Ref. (4f| ; the 
only change is in the coefficient Sj(k), which now reads 

c (lr\ _ 8 sin kl ~ sin 2fc ' iA-Ti 

6<£ 2£ k cosh(£ k /2T)/ { j 

The finiteness of the Brillouin zone results in lattice arte- 
facts in spectral functions. In particular there are cusps 
at ui = 2_Ek, when k = (71", 0,0) and (it, 7r,0). The maxi- 
mal energy is given by lu = 2Ek, when k = (tt, tt, tt). For 
a T m = 0.1, this corresponds to cusps at a T uj = 0.72 and 
1.14, and a maximal energy of a T u — 1.48. 
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